Structure formation of surfactant membranes under shear flow 
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Shear-flow-induced structure formation in surfactant-water mixtures is investigated numerically 
using a meshless-membrane model in combination with a particle-based hydrodynamics simulation 
approach for the solvent. At low shear rates, uni-lamellar vesicles and planar lamellae structures are 
formed at small and large membrane volume fractions, respectively. At high shear rates, lamellar 
states exhibit an undulation instability, leading to rolled or cylindrical membrane shapes oriented 
in the flow direction. The spatial symmetry and structure factor of this rolled state agree with 
those of intermediate states during lamellar-to-onion transition measured by time-resolved scatting 
experiments. Structural evolution in time exhibits a moderate dependence on the initial condition. 

PACS numbers: 83.50.Ax, 83.80.Qr, 87.16.dt, 87.15.A- 



I. INTRODUCTION 

Surfactant molecules in water self- assemble into vari- 
ous structures such as micelles and bilayer membranes, 
which display a rich variety of rheological properties un- 
der flowP Even if a basic structure remains to be a bi- 
layer membrane, its mesoscale structure can assume sev- 
eral different states, such as fluid Lq, or ripple P^ phase. 
Under shear flow, lamellae can be oriented parallel or 
perpendicular to the shear-gradient direction. Diat and 
Roux first discovered 20 years ago closely-packed multi- 
lamellar vesicle (MLV) structures, so-called the onion 
phase, in nonionic surfactant- water mixtures under shear 
flow.'^'Hl In the last two decades, this onion structure has 
been studied experimentally using light P^H^neutron p'^ ' ^ ' 
and X-ray scattering,i^ ^^ and also by freeze- fracture 
electron microscopy,'^ and the rheo-NM R method .'^^^^ 
Its rheology has been of large interest.^ 6 | 8 | 9 | i4 | i5 j Typi- 
cally, a critical shear rate 7c separates the lamellae and 
onion phases, where the latter phase consists of mono- 
disperse onions containing hundreds of lamellar layers. 
The onion radius ^^(7) is reversible and can be described 
by a unique decreasing function of the shear rate 7.^^ 

Time-resolving small-angle neutron scattering experi- 
ments have revealed that a two-dimensional intermediate 
structure is formed during the lamellar-to-onion transi- 
tion with increasing shear rate.^ A cylindrical or wavy 
lamellar structure was speculated to be the transient in- 
termediate structure, but could not be distinguished from 
the scattering pattern alone. Recent small- angle X-ray 
scattering experiments with increasing temperature at 
constant shear rate also indicate a similar pattern around 
the lamellar-to-onion transition. ^^ Thus, there are some 
experimental evidences of a transient state, but its struc- 
ture is still under debate. An alternative experimental 
approach to gain insight into the structural changes is 
to characterize defects observed in the lamellar stat e for 
moderate shear rates, both in surfactant membranes^^^ 



and in thermotropic liquid crystals. ^^ ^^ It is also worth 
mentioning that stable cylindrical structures on a ten- 
micrometer length scale are observed when strong shear 
flow is applied in the lamellar-sponge coexistence state.^ 

Several theoretical attempts have been made to tackle 
this complex problem of structural evolution under shear 
flow, which consider either instability of the lamel- 
lar pha se du e to undulations^^ ^^ or the break-up of 
droplet s.'^^^ Recently, a "dynamical" free energy of 
MLV under shear flow has been proposed, ^^ which takes 
into account the slow modes induced by the solvent be- 
tween the membranes together with their bending and 
stretching forces. The scaling relations for the MLV size 
and the terminal shear rate are predicted in agreement 
with the experiments of Refs. 2 3, In these theories, while 
the hydrodynamic effects of the solvent are taken into 
account, the analysis is performed for geometrically sim- 
ple structures, asp spherical onions or planar lamellae. 
Thus, the kinetic process of the transformation from the 
lamellar to the onion phase could not be investigated 
theoretically so far. In this paper, we study the detailed 
structural evolution of surfactant membranes under sim- 
ple shear flow using large-scale particle simulations. 

A few simulations have been performed for the for- 
mation of lamellar phases in shear flow, while onion 
formation has never been addressed so far. Oriented 
lamellae have been obtained in simulations of a coarse- 
grained molecular model for lipids 5^^^^^ while defect dy- 
namics has been investigated in simulations of a phase- 
field model of a smectic-A system. ^^ Onion and in- 
termediate states have large-scale structures of the or- 
der of micrometers, which are beyond typical length 
scales accessible to molecular dynamics (MD) simula- 
tions of coarse-grained surfactant molecules. In our 
study, we employ a meshless-membrane model,^^^^ 
where a membrane particle represents not a surfactant 
molecule but rather a patch of bilayer membrane - 
in order to capture the membrane dynamics on a mi- 



crometer length scale. This model is well suitable to 
study membrane dynamics accompanied by topologi- 
cal changes. Alternatively, membranes can be mod- 
eled as triangulated surfaces, ^^ 1 ^^ * ^^ * which require how- 
ever dis crete bo nd reconnections to describe topological 
changes .l^^l^^^ After the first meshless- membrane model 
was proposed in 1991,^^ sever al m eshless-membrane 
models have been developed.^^^^^^^^SU j^ contrast to 
other meshless-membrane approaches, our models^^^^ 
are capable of separately controlling the membrane bend- 
ing rigidity hc and the line tension T of membrane 
edges. Previouslyj ^^ * ^^ * we have combined our meshless- 
membran e mo del with multi-particle collision (MPC) 
dynamics,'^^^^ a particle-based hydrodynamic simulation 
technique. With MPC, the hydrodynamic interactions 
are properly taken into account, but due to the frictional 
coupling of membrane and solvent, solvent particles can 
penetrate through the membrane. We here extend the 
meshless-membrane model into explicit solvent simula- 
tion model, in which the fluid particles interact with each 
other and with the membrane via short-ranged repulsive 
potentials, so that the solvent can hardly penetrate the 
membrane, and simulate it with dissipative particle dy- 
namics (DPD), another hydrodynamics simulation tech- 
nique. 

The simulation model and methods are introduced in 
Sec.jnj Then basic membrane properties, including bend- 
ing rigidity and line tension, are described in Sec. Ill In 
Sec.llVl 



^ ^ structure formation in surfactant- water mixtures 

under shear flow is investigated for variety of shear rates 
7 and membrane volume fractions cp. At high 7 and high 
(^, a novel structure of rolled- up membranes is found, 
which are oriented in the flow direction. A summary and 
some perspectives are given in Sec. IVl 



II. MODEL AND METHODS 



All neighbor particle pairs interact via the short- 
ranged repulsive potential 



U, 
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V^ij 



in, > r-p) 
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where r^j is the distance between particles i and j. This 
potential is cut off at a distance rJ!®P = 3.2cr. The length 
a, representing the particle diameter, is employed as the 
length unit. The constant B is chosen such as to ensure 
the continuity of the potential at r = rJ^P. The (dimen- 
sionless) potential strength is set to Cc = 4. 

To favor the assembly of membrane particles into 
smoothly curved sheets in three-dimensional (3D) space, 
membrane particles interact via the additional potentials 
/7att and /7a, which have been introduced in the implicit- 
solvent version of the model previously. ^^^^ With these 
potentials, the membranes particles self-assemble into a 
single-layer sheet, which is a model representation of a bi- 
layer membrane. Here, the attractive interaction is given 

by 



f/att = 0.25 ln{l + exp[-4(pi - p*)]} - C, 



(3) 



which is a function of the local density pi of the mem- 
brane particles defined by 



Pi 



^^fcutirij/a). 
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C = 0.25 ln[l + exp(4p*)] is chosen such that /7att(0) = 
0. The cutoff function /cut in Eq. Q is a C^ function 
represented as 



/cut (5) 
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A. Coarse- Grained Model and Interaction 
Potentials 

To simulate the structure formation in surfactant- 
membranes systems, we employ a meshless-membrane 
model with explicit solvent. In this model, two types 
of particles A and B are employed, which denote mem- 
brane and solvent particles, respectively. The number of 
these particles is Nj^ and N^^ which defines the parti- 
cle density (j) = (A^^ + Nb)/V - where V is the volume 
of the simulation box - and the membrane volume frac- 
tion cp = Nj[/{Nj[ + N]s)' The particles interact via a 
potential [/, which consists of repulsive, attractive, and 
curvature interactions, t/rep, ^att, and Ua, respectively. 
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Here, the former sum for repulsive interactions is taken 
over all pairs of particles, the latter only over the mem- 
brane particles. 



where n = 12, a = ln(2){(5cut/5haif)'' 
and 5haif = 1.8 are used. We set p* 
membrane in an explicit solvent. 
The curvature potential 



Ua = ^api(ri) 



-1} with Scut = 2.1 
6 to study a fluid 



(6) 
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is introduced to incorporate the membrane bending rigid- 
ity. Here, the aplanarity api provides a measure for the 
degree of deviation of membrane particle alignment from 
a planar reference state. It is defined by 



api 



9D^ 



9A1A2A3 



T^M^ (Ai + A2 + A3)(AiA2 + A2A3 + A3Ai)' 

(7) 
where Ai < A2 < A3 are the three eigenvalues of 
the local gyration tensor aa(3 of the membrane near 
particle i, which is defined as aa(3 = ^j^j[{<^j — 
(^g)Wj - f^G)wcv{rij), with a,/3 G {x,y,z}. Here, tq = 
^jeA'^j'^cv{rij)/^j^j^Wcv{rij) is the locally weighted 



center of mass, and i^cv(^ij) is a truncated Gaussian func- 
tion 



rin 




(8) 



which is smoothly cut off at Vij = Tec The constants are 
set as n = 12, rga = 1.5cr, and Tcc = 3.2a. 

In all our simulations, the volume is chosen such that 
the number density (j) of the particles is constant, 

-3 



(j) = N/V = 0.64cr-^ N = Na^Nb 



(9) 



For higher solvent densities, the system is closer to the 
melting point, and strong attractive interaction between 
membranes sheets have been found (see Appendix A for 
details). The density (j) = 0.64cr~^ is chosen in order to 
avoid such an attraction. 
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FIG. 1: Area dependence of surface tension 7s for the ex- 
plicit solvent meshless membrane model is plotted for k^ = 
2.5,5,10, and 20 with e = 4. 



B. Thermostats 

We simulated the membranes in the NVT ensemble 
under shear flow. To keep the temperature constant, 
we employ the dissipative particle dynamics (DPD) 
thermostat,E2^^in which friction and noise forces are ap- 
plied to the relative velocities of pairs of neighboring par- 
ticles. Thus, linear and angular momentum is conserved, 
which implies that the system shows hydrodynamic be- 
havior on sufficiently large length and time scales. The 
equation of motion for the ith particle is given by 



m- 



dt 



dvi 



where Vj 



-'^{-Wij{vi-Vj)-rij^./wi]^ij{t)}rij, 

(10) 
■j/vij. Here, the weight function Wij 



conffim this dependence of surface tension, line tension, 
and bending rigidity on the control parameters for the 
explicit-solvent model. 

First, a planar fluctuating membrane is simulated with 
the particle numbers 7V^ = 1600 (and TV = 48 000 or 
64 000). For various membrane projected areas A^y = 
L^Ly, the surface tension 



% = {Pzz-{Pxx^Pyy)/2)L, 



(11) 



is investigated. Here, P„j^ is the pressure tensor given by 



N 
i=l 






'V, 



(12) 



where {/i, v} G {x, y, z}, Vi = {vf ^ v\^ v?), and the sum is 

taken over all the particles including the solvent compo- 

is Wijixi'j) = jO{Aa - Tij) with A = 2.7, where e{r) nent. In calculating P^^, the periodic image /i^ + nL^ 



is the Heaviside step function. This type of weight 
has also been used in the Lowe- Andersen thermostat.^^ 
The scheme is discretized with Shardlow-Sl splitting 
method, ^^ whose time step At^ = 0.2to is different from 
that for contributions from the molecular interactions, 
where At = O.OOSto (to = 771/7 is simulation time unit). 
In the following, we use 7 = ^Jvnk^T ja^ where k^T is 
the thermal energy unit. Although the DPD thermostat 
is usually employed for soft interaction potentials, ^^^^ it 
can also be employed for systems with steeper potentials, 
such as the Weeks- Chandler- Andersen potential. ^^ ^^ 

Shear flow with velocity v^ = 7^ in the x-direction and 
gradient in the 2:-direction is imposed by Lees-Edwards 
boundary condition. The code is optimized for use on a 
parallel computer architecture by domain decomposition 
(see Appendix B). 



III. MODEL PROPERTIES 

The meshless-membrane model with explicit solvent 
is expected to exhibit very similar equilibrium proper- 
ties as the original implicit-solvent version. ^^^^ We now 



nearest to the other interacting particles is employed, 
when the potential interaction crosses one of the peri- 
odic boundaries. 

Figure [1] shows the dependence of 7s on the projected 
membrane area A^y for ka, = 2.5, 5, 10, and 20, with 
e = 4. For 7s ^ 0, the intrinsic area A is larger than 
Axy due to the membrane undulations; buckling of the 
membrane occurs for A < A^y (the flat region at 7s < in 
Fig. [1]).^^ For tension- less membranes with Nj^ = 1600, 
the projected area A^y is given by A^y = oP^yNj, with 
aly = 1.1, 1.27, 1.33, 1.35 for k^ = 2.5, 5, 10, 20, re- 
spectively. The projected membrane area increases with 
increasing ka, because both membrane bending fluctua- 
tions and protrusions are suppressed at larger ka- 

Figure [2] shows the bending rigidity />: as a function 
of ka- Here, the bending rigidity is estimated from the 
height spectruni^ 



{\h{q)\') 



knT 



7sg^ + Kq"* 



(13) 



of the tension- less membrane (7s = 0). In calculating 
(|/i((7)P), the raw positional data of the membrane par- 
ticles (r^, z G A) is employed. ^^ ^^ Because of the slow 
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for /cc^ = 5. When e exceeds the value where F satu- 
rates — a value which becomes larger for larger ka. — 
the membrane particles prefer to reside at the edges, be- 
cause the curvature force is not strong enough to avoid 
aggregation due to the stronger attractive forces. 

For our simulations of membrane ensembles in shear 
flow, we choose the model parameters /Cc^ = 5 and e = 4, 
where the membrane has bending rigidity n/kBT = 
11 ± 1. With the estimated values for the line tension 
Tcr /k^T = 4.2 ± 0.2, the relaxation time scale of struc- 
tural transitions can be characterized by 
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20 
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(15) 



FIG. 2: Dependence of the bending rigidity k on ko, is plotted 
for e = 4, estimated with the use of a t ension-less planar 
membrane at Na — 1600 through Eq. (13). As shown in the 
inset, the height spectrum of the membrane exhibits a ^~^ 
spectrum. 



where 77 is the solvent viscosity, and Re is the critical 
radius of a flat disk. Assuming a flat disk with radius R 
and a vesicle with the same membrane area, we obtain 
the corresponding elastic free energies J^d = 27tRT and 
Ts = 87r(A>: + R/2) c^ Attk.; thus, a disk exhibits transition 
to closed vesicle shape if the radius is around 




Rc = 2i^/r ^ (5.3±0.5)cr. 



(16) 



FIG. 3: e dependence of the line tension F of explicit-solvent 
meshless-membrane model, calculated for /c^ = 5 and 10 with 



the use of Eq. ( 14 ) 



dynamics of long-wavelength height fluctuations, it be- 
comes more time-consuming to obtain precise data than 
for the implicit-solvent model. ^^ ^^ Therefore, {\h{q)\'^) is 
measured using 16 independent runs for each /c^, and the 
averaged spectrum is then fitted to Eq. ([l3| in the range 
q < 1.2a~^. As in the implicit model, the bending rigid- 
ity is found to be proportional to k^^ which demonstrates 
the controllability of hi in our model. 

In Fig. [3) the e dependence of the line tension F is 
displayed for /cq, = 5 and 10 for a membrane strip with 
two edges with lengths equal to L^ with A^^ = 1600. 
Here, F is determined via the relatioiP^ 



F = {{Pyy + P,,)/2 - P,,)LyLj2. 



(14) 



At a small e, F is proportional to e and almost inde- 
pendent of /cq,, similarly to the implicit-solvent meshless 
membrane model^^. While F increases linearly up to 
e = 8 for /ccK = 10, it levels off and saturates at e = 4 



With the solvent viscosity rj = (2.1 ± 0.1) x m{(jto)~^ 
obtained from a solvent-only simulation, an estimation 
of the relaxation time yields r = 28.4to. In the follow- 
ing, the time will be measured in unit of r. Since the 
membrane thickness, typically around 5nm for non- ionic 
surfactants like polyethylene-glycol-ethers CnE^, corre- 
sponds to the size a of the membrane particles in our 
simulation, r is equivalent to about 0.36/is (with the vis- 
cosity of water at 300K, r]^ :^ 0.8mPa • s). 



IV. STRUCTURE FORMATION WITH AND 
WITHOUT SHEAR FLOW 



We now employ the meshless-membrane model 
with explicit solvent to study structure formation in 
surfactant- water mixtures, both with and without shear 
flow. In all simulations, the total particle number is fixed 
at A/" = Nj[ + Njs = 960 000, and thus, the system size 
is a cubic box with side length L = 114.47cr. Simula- 
tions are performed for various membrane volume frac- 
tion cp = Na/N, with cp = 0.0625, 0.125, 0.1875, 0.25, 
and 0.3125. The dynamical evolution is integrated over a 
total time interval of 1.2 x 10^ MD steps, corresponding to 
2.11 X lO^r. All the particles of both species are initially 
distributed randomly in the simulation box. Averages 
are calculated over the last 2 x 10^ steps, where the sys- 
tem is assumed to have reached a stationary state. After 
briefly explaining the structures obtained by equilibrium 
simulations without shear in Sec. |IVA] we present results 
for the structure formation in a system under linear shear 
flow in Sees. HVBl and HVCl 



(a) 7 = 0,93 = 0.0625 



(b) 7 = 2.84 X 10"^ (^ = 0.0625 (c) ^r = 5.68 x 10^^ ip = 0.125 




(d) 7r = 2.84 x 10"^ 9? = 0.1875 (e) 7r = 1.42 x 10-\ if = 0.1875 (f) 7r = 5.68 x 10-\ ip = 0.1875 






FIG. 4: Snapshots of the configuration of membrane particles for cp = 0.0625 (a) without shear (7 = 0) and (b) under shear 
flow with shear rate jr = 0.0284, and (c) for cp = 0.125 with jr = 0.0568. Snapshots for cp = 0.1875 under shear flow are 
shown with shear rates (d) jr = 0.0284, (e) jr = 0.142, and (f) jr = 0.568. In (b)-(f), the (average) imposed flow velocity is 
V = jzex. Solvent particles are not displayed. 



A. Mesophases in Thermal Equilibrium 



At a low membrane volume fraction cp = 0.0625, 
membrane particles self- assemble into vesicles, each of 
which is composed of around 100 membrane particles (see 
Fig. |4ja)). This result is consistent with Eq. (16), which 
predicts the critical particle number of a vesicle to be 
Nc = TrR^/a^ c^ 75. At a higher membrane volume frac- 
tion (f = 0.125 {Nji = 120 000), the membrane surfaces 
is found to percolate through the whole system via the 
periodic boundaries. This behavior is not unexpected, 
because each disk with radius Re covers a region of vol- 
ume Vc = {2Rc)^ by rotational diffusion, and therefore 
these disks should overlap and merge when more than 
^c ^ y/'^c = 1200 vesicles are present, which exceeds the 
number of vesicles with critical size (Nj^/Nc). 



of the membrane density is calculated as 



For cp = 0.25 and 0.3125, periodic lamellar states are 
formed owing to the repulsive interactions between the 
membranes, as shown in Fig. Isja). The membranes are 
curved to fill up space with random orientations. In 
thermal equilibrium, these lamellar layers would prob- 
ably have a unique orientation throughout the system; 
this well-ordered state is difficult to reach in simulations 
because the structural relaxation time well exceeds the 
accessible simulation time scale. The 3D structure factor 
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dr 
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where Shj{{r) = ^i^j^cr^^{r — Vi) — (j) is the local devi- 
ation of the membrane particle density from its average. 
Figure |5|b) shows Sj^{q) for Lp = 0.3125. The scatter- 
ing intensity is spherically symmetric, as demonstrated 
by a two-dimensional (2D) color-map for q^ = 0. Peaks 
arising from the inter-lamellar distance are observed at 
\q\ = qi = 1.74cr"^ and (72 = 3.49cr"^ = 2qi with heights 
9.8 and 0.78, respectively. The former corresponds to 
length of L = 2i\ jqx = 3.61cr, which provides a precise 
estimate of the interlayer distance. 



B. Dynamic Phase Diagram under Shear Flow 

In shear flow, the vesicle and lamellar states depend on 
the concentration cp and the shear rate 7, as displayed in 
Fig. |6] At (p = 0.0625, assemblies of uni-lamellar vesicles 
are observed for 7r < 0.02. At a higher shear rate, the 
membrane particles tend to assemble into plate-like mem- 
brane disks, which then align parallel to the shear flow 
direction, as demonstrated by the comparison of snap- 
shots in Figs. [4] (a) and (b). At (f = 0.125, lamellar lay- 
ers with non-uniform lamellar distances are observed for 






% 



JT < 0.06, as shown in Fig.Qc). There is a larger variety 
of phases at (^ = 0.1875; at low shear rates jr < 0.1, the 
system is a mixture of lamellae and cylinders as shown in 
Fig.Qd); at jr = 0.142, the lamellae are rolled up collec- 
tively, see Fig. [4]^e); finally, at large shear rates jr > 0.3, 
the system exhibits a reentrant lamellar state, as shown 
inFig.gf). 

At higher cp, the lamellar states align perpendicularly 
to the shear-gradient (z) direction in the regime of low 
shear rates {jr < 0.1). At larger 7, they exhibit an in- 
stability to a rolled-up shape whose axis is parallel to the 
flow direction, which will be investigated in more detail in 
Sec. |IVC| below. At cp = 0.25, there is again a reentrant 
behavior of lamellar states, i.e. nearly planar aligned 
layers appear at large shear rates jr > 0.284. Note 
that experimental phase diagrams often exhibit reen- 
trant behaviors^ ^ with increase in the shear rate at a 
certain membrane volume fraction; the mixture changes 
from lamellar state to onion state, and then after going 
through the coexistence region it enters again into an 
oriented lamellar state again. Although the onion state 
with densely packed MLVs has not been obtained in our 
simulations, the reentrant behavior is qualitatively con- 
sistency with the experimental observations. 



C. Rolled-up Lamellar Structures 



FIG. 5: (a) Snapshot of the configuration of membrane par- 
ticles for (p — 0.3125 without shear fiow (7 = 0). To facili- 
tate visualization, only a thin planar slice is displayed. Sol- 
vent particles are not shown, (b) Structure factor Sa[q = 
(0,^^,^^)] for (^ = 0.3125 at 7 = 0. 
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FIG. 6: Dynamic phase diagram of the explicit-solvent 
meshless-membrane model as a function of the volume frac- 
tion cp of the membrane component and the shear rate 7. The 
parameters are /c^ = 5, e = 4, = 0.64cr~^, and N — 960, 000. 
The dashed lines are guides to the eye. 



1. Structure Analysis 

Snapshots of membrane conformations are shown in 
Fig. [7] for (p = 0.3125 at various shear rates jr = 0.0284, 
0.0568, and 0.142. In ah simulations for (p > 0.125, the 
membranes are completely aligned with the fiow direc- 
tion. Thus, as shown in the bottom panel of Fig. [7| 
the membrane configurations can be visualized by cross- 
sectional slices perpendicular to the fiow direction. Fig- 
ure [Jl demonstrates the transition from the lamellar state 
to the rolled-up state, which is stable in the region 
(p > 0.175 and 7r > 0.1 in the phase diagram of Fig. |6] 

The structural changes accompanying this instability 
can be characterized by the average orientation and the 
mean-square local curvature of the membrane. The nor- 
mal unit vector n is calculated fr om the first-order mov- 
ing least-squares (MLS) methocP^^^^ applied to the 
configurations of membrane particles. Using a weighted 
gyration tensor aa/s = Y^ji^'j - ^G')(/^j " /^G)^cv(nj), 
where a, /3 G {x^y^z}^ n is obtained as an eigenvector 
corresponding to the minimum eigenvalue of aQ,/3, which 
together with the other two eigenvectors ei and €2 con- 
stitutes an orthonormal basis. Here, the cut-off function 
'^cv(^ij) in Eq. (I8| is employed with the same cut-off 
lengths r^^ = 3.2(7. As a quantitative measure for the 
undulation instability, we employ the orientational order 
parameter 



S, = [2{n • e,f - 1] = cos(2^) 



(18) 



(a) 7r = 2.84 x 10" 



(b) JT = 5.68 X 10- 





(c) 7r = 14.2 X 10-2 







FIG. 7: Snapshots membrane conformations for volume fraction cp = 0.3125, with shear rates (a) jr — 0.0284, (b) jr = 0.0568, 
and (c) JT = 0.142. Views from the flow (x) direction are shown in the upper panels. Corresponding cross-sectional views 
(with —3.0(7 < X < 3. Oct) are shown in the lower panels. 
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FIG. 8: (a) Membrane orientational order parameter (Sz) in 
Eq. ([18| for (^ = 0.125, 0.1875, 0.25, and 0.3125, and (b) the 
mean-square local curvature {H'^) for cp = 0.1875, 0.25, and 
0.3125, both as function of the shear rate jr. 



with the normalization for the two-dimensional order of 
cylindrical and planar symmetry. 

The instability can also be characterized by calculat- 
ing the mean-square curvature. The second-order MLS 



method^^ provides an estimate of the membrane curva- 
ture from the particle configurations in the following way. 
For each particle i, we perform a rotational transforma- 
tion into the principal coordinate system of the gyration 
tensor of the neighbor particles j around i's weighted 
center of mass re by 







rof, 



and then employ the parabohc fit function 
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where the coefficients of the Taylor expansion 
2^0, hx^ hy, hxx^ hyy aud hxy are fitting parameters.^ 
By a least-squares fit, the estimated value of the mean 
curvature H = (Ci +C2)/2 for particle i is then obtained 
as 
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Figure [S] displays the results for the spatial average 
{Sz) of the 2D orientational order parameter and for 
the mean-square local curvature (i^^). When the mem- 
branes are aligned with the x — y plane orthogonal to 
the shear-gradient direction, Sz becomes unity. As the 
membranes roll up, Sz decreases. Here, (Sz) = for 
perfectly cylindrical state (where (n • e^^)^ = and 



4 
2 

qz 

-2 

-4 



(a) 7r = 2.84 x 10" 



'iii-if' 



(b) 7T = 5.68 X 10" 



(c) 7r = 14.2 X 10" 



-4 









^1/ 







n200 
100 




FIG. 9: Color maps of the structure factor S a{S^ ^ Qy ^ Qz) for volume fraction Lp — 0.3125 and shear rates (a) 7T = 0.0284, (b) 
7T = 0.0568, and (c) 7T = 0.142. Since the structure is uniform in the flow (x) direction, only data for g^ = are shown. 



{n-CyY = {n-ezY = 1/2 ), and (Sz) = —1 for a perfectly 
flat lamellar layers perpendicular to the vorticity {y) di- 
rection (where (n • e^)^ = 0). For all (p > 0.1875, the 
values of {Sz) and {H'^) provide evidence for rolled- up 
structures at jr = 0.142. For (p = 0.1875 and 0.25, (Sz) 
increases again at ^yr = 0.284, indicating reentrant align- 
ment of the membranes in a lamellar stack. However, for 
(f = 0.3125, (Sz) remains small (and {H'^) large), which 
implies that rolled- up structures exist also at jr = 0.284. 
Thus, the evolution of rolled- up conformations is the 
most pronounced at the highest membrane volume frac- 
tion (f = 0.3125. 

In Fig.|9J results for the structure factor Sj[{q) at q^ = 
are shown for the same set of data as in Fig. [7| Due to 
the nearly complete alignment of membranes in the flow 
{x) direction, all the structural features are reflected in 
SA{q) in the Qy — Qz wave- vector plane. Since lamellar 
layers are nearly planar at the low shear rate ^yr = 0.0284, 
a sharp peak is observed around {qy^ Qz) = (27rg'i, 0). Be- 
cause of the undulation instability, the pattern becomes 
more circular at higher 7. 

In the small-angle neutron^ and X-isiy^^ scattering ex- 
periments, the scattering beams can be injected from two 
directions, either radial or tangential to the shear cell, 
which correspond to shear-gradient and flow directions, 
respectively. After the shear flow is applied, after a while 
a Bragg peak of the radial beam develops in the vorticity 
direction, while the scattering pattern becomes isotropic 
for the tangential beam. This suggests 2D-isotropic un- 
dulations of the lamellar structure perpendicular to the 
flow direction. Later in time, the radial beam is scattered 
isotropically in the tangential direction, indicating the 
formation of onion structures. The structure factor S{q) 
in our simulation (Fig. |9|c)), agrees well with S{q) of 
this transient states in these scatting experiments. Mea- 
surements of solvent diffusion in a Rheo-NMR experi- 
ment of a non-ionic surfactant systenP^ show a diffu- 
sion anisotropy in the direction of shear flow in the inter- 
mediate state. These experimental results suggest that 



the membranes are aligned in the direction of shear flow 
in the intermediate state of the lamellae-to-onion tran- 
sition, although more detailed information on the struc- 
tural membrane arrangement could not be obtained. The 
rolled-up lamellar structures observed in our simulations 
match the experimental evidence, so that they are a good 
candidate for the intermediate states. 



2. Temporal Evolution of Membrane Structures 

In simulations of a large system, even if extended over 
a long time interval, the resultant structure often ex- 
hibits dependence on the initial conditions - owing to the 
slow processes involved in the dynamics. Therefore, we 
compare here the time evolution at (p = 0.3125 starting 
from both, a lamellar state and a random distribution 
of membrane particles (the latter corresponding to the 
simulations described in the previous subsections). For 
an initial lamellar state, the final configuration of a sim- 
ulation run with 1.2 x 10^ MD steps at a small shear rate 
jr = 0.0284 (see Fig. JTJ^a)), is employed. 

In Fig.fToJa), the orientation order parameter (Sz) for 
both initial conditions are shown for jr = 0.568, 0.284, 
and 0.142. For the case of a random initial distribu- 
tion, small disks merge into randomly oriented surfaces, 
which then align in the shear flow to become a nearly 
perfect lamellar stack with some defects. Afterwards, 
lamellae roll up into slightly larger rolls. During rolling 
up? (Sz) exhibits an overshoot (see Fig. 10 ^c)), and fi- 
nally approaches a constant as the structure relaxes into 
a (meta) stable state. The overshoot amplitude depends 
on initial states and random noises. 

Figs. floFb) and (c) show the membrane conformations 
after an elapsed time of t = 2.11 x lO^r (for jr = 0.284) 
for the two types of initial conditions, and explain the 
origin of the substantially different values of the order 
parameter (Sz) in Fig. fToJa). When the random state is 
taken as initial condition, rolled-up structures are con- 
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FIG. 10: Comparison of structure formation from a random 
initial configuration and from a lamellar state at (/? = 0.3125. 
(a) Time evolution of the orientational order parameter {Sz)- 
The dotted and solid lines show the data starting from the 
random and lamellar initial states, respectively. Blue, red, 
and black lines represent the shear rates ^r — 0.568, 0.284, 
and 0.142, respectively. Snapshots of the final configurations 
for 7T = 0.284 at t = 2.11 x lO^r, as they have developed 
from (b) the random and (c) the lamellar initial states, are 
also shown. 



siderably more pronounced; this may be traced back to 
the presence of defects in the lamellar structure which 
forms at short times tjr ~ 100. 



V. SUMMARY 

In this paper, we have constructed an explicit-solvent 
meshless-membrane model for surfactant- water mixtures. 
The model reproduces properties of an earlier implicit- 
solvent meshless membrane model, where membrane 
bending rigidity and line tension can be independently 
controlled to a large extent. The model enables large- 
scale simulations of structural changes, where dynamical 
effects of hydrodynamic interactions have to be taken 
into account. At present, such a large simulation with as 
many as one-million particles can be realized by paral- 
lelized molecular dynamics simulation methods. 

Our main results concern the effects of shear flow on 
the structure formation of membrane ensembles. Various 
structures including vesicles, lamellae, and multi-lamellar 
states with nearly cylindrical symmetry have been found, 
most of which are qualitatively consistent with experi- 
mental observations of non-ionic surfactant membranes 
under shear flow. Especially, a cylindrical instability of 
multi-lamellar membrane is predicted to occur perpen- 
dicularly to the flow direction. The corresponding scat- 
tering patterns are in qualitatively agreement with the 
results of small angle neutron (and X-ray) scattering ex- 
periments under shear deformation. The rolled-up lamel- 
lae are a good candidate for the intermediate structures 
on the way to the onion state, which are observed in the 
experiments. 

Our simulations do not reproduce onion formation, 
which is ubiquitously observed in experiments on /im 
length scales. We speculate that it is due to the lim- 
ited system size, which can be overcome by larger-scale 
calculations in the future. On the other hand, in the 
experiments, strains larger than ^t > 10^ are necessary 
to reach the onion state, which indicates that very long 
simulation runs are required to obtain these states. 

The control of the physical parameters including the 
line tension T, the bending rigidity k^ and the Gaus- 
sian modulus K is another challenge. While T and k are 
easily controlled in our model, k is more difficult. The 
Gaussian modulus K might also play an important role 
in the structure formation, because it is directly related 
to topological changes happening on the way of onion 
formation. ^^ 



For the case of a lamellar initial configuration, the un- 
dulation instability becomes more conspicuous when the 
applied shear flow is stronger. Moreover, while strong 
undulations are observed at low 7 with the random ini- 
tial configuration, less undulations take place with the 
lamellar initial configuration, as indicated by {S^) ^ 1. 
Thus, the initial conditions play an important role in the 
selection of the transition path and structure formation. 
It may depend not only on the shear rate and relaxation 
time but also on the distribution of structural defects in 
the lamellar layers. Thus, more systematic studies are 
required in the future to clarify the hysteresis of these 
systems. 



Appendix A: Solvent-Mediated Forces at Higher 
Solvent Density 

In the present simulation model, solvent particles have 
a similar size as membrane particles. Here, we discuss 
the finite-size effects of the solvent particles, if much 
higher solvent densities are used than in the present 
study. When the solvent particles are densely packed, 
the system approaches a crystallization transition, and 
local crystalline order emerges to bring about interac- 
tions between closely spaced lamellar layers. 

As an example, we study the explicit-solvent meshless- 
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FIG. 11: Snapshot of the MLV state due to the solvent- 
mediated attractive interactions for the system <S2 at = 
0.72, LP = 0.06 and N = 480,000 with different size ratio 
(^bI c^A — 1-2, for shear rate 7r = 0.0142. 



FIG. 12: Equihbrium lamellar distance diam versus initial 
lamellar distance dini between a tension-less planar membrane 
and two (smaller) membrane discs. For the systems S\ and 
<S2, lengths are represented in units of a and g^^ respectively. 
The dotted line indicates the line dini = c/iam- 



membrane model for higher number density (denoted sys- 
tem (52), where (\) = {(t\Na^(tINb)/V = 0.72 with the 
total particle number N = N^ ^ Nb = 480 000). The 
particle radii of the two components are chosen such as 
as = 1.2(7^, where a a and as denote the radii of mem- 
brane and solvent components, respectively. The cut-off 
lengths of the interactions are set to rJ!®P = 2.7<j^, rga = 
1.5(7^, Tec = 3.0(7^, respectively. The repulsive in- 
verse twelfth-power potential exhibits melting transition 
around the volume fraction around 0.43 (corresponding 
to (j) 0:^ 0.8 in our definition). ^^ Thus, our density ^ = 0.72 
is close to the crystallization line. 

As an example. Fig. [TT] shows a snapshot for membrane 
volume fraction cp = Nj^a\/{Nj,a\ + Nscr^) = 0.06 for 
a shear rate ^yr = 0.0142. After an initial relaxation 
from a random configuration of the membrane particles, 
the membrane layers assemble into stacks, and then form 
MLVs after gathering a certain amount of the membrane 
patches. In the snapshot, the layers are stacked with dis- 
tances d c^ 3.0(723, which seems to arise from the discrete- 
ness of the interstitial solvent particles - here it should 
be noted that interlayer distance cannot be smaller than 
d 0:^ 2.0(723 because it is inside the range of Urep and /7q,, 
which both act repulsively. 

To avoid this solvent-mediated force, a longer cut-off 
range of the repulsive forces and a lower density are em- 
ployed in the simulation described in the main text (de- 
noted system (51)). To demonstrate the difference be- 
tween the systems 51 and 52, we simulate three lamellar 
layers in the following way. A periodic simulation box 

with lengths L^ = Ly = ^J^y/(JA, ^z = V/{L^Ly) 
is set up so that a tension-less membrane (with particle 
number A^^o) = 1000 is put along the x?/-plane. Here, 
the projected areas are set to A^^y = 1.267A^^o for 51 and 
A^y = 1.33bNAo for 52, respectively. With various inter- 
layer distance dini, circular membrane disks with particle 
number A^^i = Na2 = 400 are put on both sides of the 



membrane. Solvent particles are then inserted to fill up 
the system; they are relaxed for fixed membrane configu- 
ration for 10^ MD steps. Afterward, a simulation of the 
full system is performed for 1.5 x 10^ MD steps to ob- 
tain an equilibrium interlayer distance (iiam- As shown in 
Fig. [T2J while the interlayer distance increases with time 
in the 51 case, there are stable lamellar distances around 
diam = 3.5(7^ and 4.7(7^ in the 52 case with higher sol- 
vent density, which correspond to 3(723 and 4(723, respec- 
tively. Thus, in 52, the effective attractive potentials lead 
to the stable multi-lamellar layers due to the attractive 
interactions. 



Appendix B: Numerics and Parallelization 

Numerical simulations have been carried out on mas- 
sively parallel supercomputers. On 256 CPUs of Intel 
Xeon X5570 (2.93GHz) in SGI Altix ICE 8400EX at 
ISSP, it costs 72 hours to perform a run of 1.2 x 10^ 
simulation steps for (f = 0.3125 and N = 960, 000. Here, 
about 10% of the total time is for the force calculation 
of t/rep, 30% for Ua, arouud 20% for the construction of 
the buffer, and 20% for the communication between the 
MPI processes. The system is divided into cubic (or rect- 
angular) boxes, each of which is calculated by one MPI 
process. We also parallelize calculations of each process 
with the use of OpenMP by performing calculations of 
different far away pairs at once on different threads. The 
code is optimized to achieve 15% performance compared 
with the theoretical limit on X5570 processors. 

Each process has separate cell lists and neighbor lists 
for solvent and membrane particles. To apply shear, we 
employ a simultaneous affine deformation of the total sys- 
tem, each MPI process box, and each cell for neighbor 
search, so that a square is transformed into a parallel- 
ogram shape consistent with the shear deformation, as 
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FIG. 13: Schematic picture in two dimensions for the imple- 
mentation of Lees- Edwards boundary conditions in a program 
parallehzed with MPI communication. 



parallelograms are reflected to make the strain -0.5, see 
Fig. [iSJ which basically requires an all-to-all communi- 
cation of all the position and velocity data. 
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